library(AER)
library(gmm)
library(MASS)
library(stargazer)

#1 preparations
dt = read.csv (file = "boston.csv",header = TRUE, sep = "," , row.names = NULL,stringsAsFactors = FALSE)
dt = dt [ complete.cases( dt ) , ]
#complete.cases return true if one row doesn't contain NA
n = nrow ( dt )
head(dt[,c( "value" , "crime" , "industrial" , "distance" , "black" , "ptratio" ) ] )

#extract variables and iv
name_x=c("constant","crime","industrial","distance")
name_z=c("constant","black" , "ptratio" ,"industrial","distance")
constant=rep(1,n)
X=cbind(constant,dt[,c("crime","industrial","distance")])
colnames(X)=name_x
Z=cbind(constant,dt[,c("black" , "ptratio" ,"industrial","distance")])
colnames(Z)=name_z
Y=dt[,"value"]

#2 step gmm
#2.1recall ole estimates for later comparison
model=value~crime+industrial+distance
ols=lm(model, data=dt)
#summary of ols
coeftest(ols,vcov.=NeweyWest)
cov_hac=NeweyWest(ols)
se_hac=sqrt(diag(cov_hac))
t_hac=coef(ols)/se_hac
p_hac=pt(abs(t_hac),df=nrow(dt)-ncol(X),lower.tail = FALSE)*2 #make p value more conservative

#2.2perform 2sgmm
Z=as.matrix(Z)
X=as.matrix(X)

ZX=t(Z) %*% X/n # or ZX = crossprod(Z,X)/n
ZY = t(Z) %*% Y/n
var_g_1 = t(Z) %*% Z/n
theta_1= solve(t(ZX) %*% solve(var_g_1) %*% ZX) %*%t(ZX) %*% solve(var_g_1)%*% ZY
cov_1=solve(t(ZX)%*%solve(var_g_1)%*% ZX)
#calculate the estimates and check
se_1=sqrt(diag(cov_1))
t_1=theta_1/se_1
p_1=pt(abs(t_1),df=nrow(dt)-ncol(X),lower.tail = FALSE)*2 #make p value more conservative
print("1st step estimates:"); theta_1
print("1st step se estimates:"); se_1
print("t statistics:"); t_1
print("p values:"); p_1
##2nd step
Zu_1=t(Z)%*%(Y-X%*% theta_1)
var_g_2=Zu_1%*% t(Zu_1)/n
u_1=Y-X%*% theta_1
i=1
l=1
L=10 #for simplicity <=n
weight=rep(0,L)
Q=matrix(0,nrow = ncol(X),ncol = ncol(X))
Q_1=matrix(0,nrow = ncol(X),ncol = ncol(X))
Q_2=matrix(0,nrow = ncol(X),ncol = ncol(X))
for (l in 1:L) {
  weight[l]=1-l/(L+1)
}
for (i in 1:n) {
  Q_1=Q_1+u_1[i]^2*(X[i,]%*% t(X[i,]))
  
}
Q_1 = Q_1/n

for (l in 1:L) {
  for (i in ( l+1) : n ) {
    Q_2=Q_2+weight[l]*u_1[i]*u_1[i-l]*(X[i,]%*% t(X[i-l,])+X[i-l,]%*% t(X[i,]))
  }
  
}
Q_2=Q_2/n
Q=Q_1 + Q_2

XXT=X%*% t(X)
omega=ginv(XXT)%*% (X%*%Q%*%t(X) )%*% ginv (XXT)
var_g_2=t(Z)%*% omega%*% Z/n
theta_2=ginv(t(ZX)%*% ginv(var_g_2)%*% ZX)%*%t(ZX)%*% ginv ( var_g_2) %*% ZY
cov_2=ginv(t(ZX)%*% ginv(var_g_2)%*% ZX)
#calculate the estimates and check
se_2=sqrt(diag(cov_2))
se_2=as.vector(se_2)
t_2=theta_2/se_2
p_2=pt(abs(t_2),df=nrow(dt)-ncol(X),lower.tail = FALSE)*2 #make p value more conservative
print("2nd step estimates:"); theta_2
print("2nd step se estimates:"); se_2
print("t statistics:"); t_2
print("p values:"); p_2

#finally theta_1 is not optimal
se_1>=se_2
#display the results 
stargazer(ols,ols,ols,coef = list(coef(ols),theta_1,theta_2),se=list(se_hac,se_1,se_2),
          t=list(t_hac,t_1,t_2), p=list(p_hac,p_1,p_2),
          t.auto = FALSE,p.auto = FALSE,
          report = "vcstp",keep.stat = c("n"),
          add.lines = list(c("Type", "HA Roubusted OLS" , "1st step GMM", "2nd step GMM")),
          title="Results of rols and numerical 2 step GMM",
          label = "Numeric",
          type = "latex",header = FALSE,font.size = "small",#generating latex code
          table.placement = "htb",omit.table.layout = "n")
